require(ggplot2)
require(ggthemes)
require(tcltk)


rm(list=ls())


#This script reads in all of the digitized data for performance metrics, then analyzes it
  #It measures performance metrics over 2 different time periods
  #In phase 1 we measure heading (yaw) change, duration, and avg. angular velocity.
    #This is the phase in which they are turning away from the stimulus, so we mostly care about how quickly they can reorient themselves
  #In phase 2 we measure COM displacement. Since we record points over a standard duration, this is colinear with COM linear velocity.
    #This is the phase in which they are escaping, so we record their ability to move their body away from the stimulus.


dir <- "Digitized midline points"
setwd(dir)


#Read the names of the .csv files in the current folder
fileNames <- list.files(pattern="*.csv", all.files=TRUE, no..=TRUE)
#Open all the files
list2env(
  lapply(setNames(fileNames, make.names(gsub("*.csv$", "", fileNames))), 
         read.csv), envir = .GlobalEnv)

#Preallocate a data frame to hold all the curvature data
finalPerfData <- data.frame(Treatment=character(),Animal=integer(),Sequence=integer(),ID=character(),yawDeg=double(),P1Duration=double(),angVel=double(),comVelMMms=double())
#This loop open each file, performs the analysis, then saves out the data in a new folder

#now run the actual analysis
for(i in 1:length(fileNames)){
  nameNow <- fileNames[i]
  #Get rid of ".csv" in the filename
  nameNow2 <- substr(nameNow,1,nchar(nameNow)-4)
  #Open the file
  Data <- eval(parse(text= nameNow2))
  
  #rename variables in data frame
  names(Data) <- c('x1','y1','x2','y2','x3','y3','x4','y4')
      #x1 and y1 are the snout
      #x2 and y2 are the back of the head
      #x3 and y3 are the COM
      #x4 and y4 are the peduncle
  
  
  #Extract video data from filename
  splitStrings <- strsplit(nameNow2,fixed=FALSE, split = '_')
  treatNow <- splitStrings[[1]][3]
  aNow <- as.integer(substr(splitStrings[[1]][4],2,3))
  TANow <- paste(treatNow,aNow,sep='')
  sNow <- as.integer(substr(splitStrings[[1]][5],2,3))
  
  
  #Make time variable (all videos are 1k hz)
  Data$time <- seq(from=1, to=nrow(Data), by=1)
  
  realData <- Data[complete.cases(Data),]
  
  
  
  #Phase 1 measurements
      
          #Heading Change
      snoutXt1 <- realData$x1[1]
      snoutYt1 <- realData$y1[1]
      snoutXt2 <- realData$x1[2]
      snoutYt2 <- realData$y1[2]
      branchXt1 <- realData$x2[1]
      branchYt1 <- realData$y2[1]
      branchXt2 <- realData$x2[2]
      branchYt2 <- realData$y2[2]
      
      deltaXT1 <- branchXt1-snoutXt1
      deltaYT1 <- branchYt1-snoutYt1
      deltaXT2 <- branchXt2-snoutXt2
      deltaYT2 <- branchYt2-snoutYt2
      
      vectT1 <- c(deltaXT1, deltaYT1)
      vectT2 <- c(deltaXT2, deltaYT2)
      
      magVT1 <- sqrt(deltaXT1^2+deltaYT1^2)
      magVT2 <- sqrt(deltaXT2^2+deltaYT2^2)
      
                      #The  dot product of two vectors is x*y = |x||y|cos(theta)
                        #So the angle between the vectors is arccos(x*y/|x||y|)
                        #The r function for dot product is %*%
      
      yawRad <- acos((vectT1%*%vectT2)/(magVT1*magVT2))
      yawDeg <- as.vector(yawRad*180/pi)
      
      
      
          #Duration of phase 1
      P1Duration <- realData$time[2]-realData$time[1]
      
          #Average angular velocity over phase 1
      angVel <- as.vector(yawDeg/P1Duration)  #This is in degrees per millisecond
  
  #Phase 2 measurements
  
      # COM Displacement
  comXT2 <- realData$x3[2]
  comXT3 <- realData$x3[3]
  comYT2 <- realData$y3[2]
  comYT3 <- realData$y3[3]

  comDispPix <- sqrt((comXT3-comXT2)^2+(comYT3-comYT2)^2)
  #now convert to cm. Calibration value is 7.5 pixels/mm
  comDispMM <- comDispPix/7.5
  
  #convert to linear velocity by dividing by 8ms
  comVelMMms <- comDispMM/8
  
  
  
  #Save all the data from this video
  finalPerfData[i,] <- data.frame(Treatment=treatNow,Animal=aNow,Sequence=sNow,ID=TANow,yawDeg=yawDeg,P1Duration=P1Duration,angVel=angVel,comVelMMms=comVelMMms)
  
}

ggplot(finalPerfData, aes(x=Treatment,y=angVel, colour=Treatment)) + 
  geom_boxplot() +
  geom_point(aes(),size=2,alpha=1) +
  theme_few()


#show max values for each individual so I can enter it into spreadsheet

allIDs <- unique(finalPerfData$ID)
for (k in 1:length(allIDs)) {
  
  temp <- max(finalPerfData$angVel[finalPerfData$ID==allIDs[k]])
  temp2 <- paste(allIDs[k],' ',temp)
  print(temp2)
}


#Put max values in data frame
    #make new columns to hold only max values
naVect <- rep(NA,nrow(finalPerfData))
naDF <- data.frame(maxYaw=naVect, maxDuration=naVect, maxAngVel=naVect, maxLinVel=naVect)
DataWithMax <- cbind(finalPerfData,naDF)

  #Find the maxima, find where they are in the data frame, and copy the value to a new column for that metric
        #doesn't work for duration because it is copying all instances that match each maximum duration value.
        #I don't care because duration isn't a performance metric, more of an effort metric.
for (k in 1:length(allIDs)) {
  for (j in 5:8){
    tempPerf <- max(DataWithMax[DataWithMax$ID==allIDs[k],j])
    DataWithMax[DataWithMax[,j]==tempPerf,j+4] <- tempPerf
    
  }
}


#I am saving this data to "/performance data/PerformanceData_DATE.csv"
fileSaveTo1 <- tclvalue(tkgetOpenFile(default=dir))
write.csv(finalPerfData, file=fileSaveTo1)

#I am saving this data to "/performance data/PerformanceDataWithMax_DATE.csv"
fileSaveTo2 <- tclvalue(tkgetOpenFile(default=dir))
write.csv(DataWithMax, file=fileSaveTo2)
